Run-and-tumble particles with hydrodynamics: 
sedimentation, trapping and upstream swimming 
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We simulate by lattice Boltzmann the nonequilibrium steady states of run-and-tumble particles 
(inspired by a minimal model of bacteria), interacting by far-Held hydrodynamics, subject to con- 
finement. Under gravity, hydrodynamic interactions barely perturb the steady state found without 
them, but for particles in a harmonic trap such a state is quite changed if the run length is larger 
than the confinement length: a self-assembled pump is formed. Particles likewise confined in a 
narrow channel show a generic upstream flux in Poiseuille flow: chiral swimming is not required. 
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The motility of microorganisms raises basic physics 
questions that range from local swimming mechanisms 
[TH5] to many-body emergent phenomena [H 03 [TJ. In 
the latter context, even grossly simplified models repre- 
sent a challenging and active area of nonequilibrium sta- 
tistical mechanics [1]. In some cases experimental near- 
counterparts to these models can be devised in which 
various complicating factors (cell division, chemotaxis, 
etc.) are environmentally or genetically suppressed [5]. 

Indeed certain bacteria, including E. coli, exhibit mo- 
tion which can be idealized as a 'run-and-tumble' model. 
Here straight 'runs' at constant speed v are punctuated 
by sudden, rapid and complete randomizations in direc- 
tion, or 'tumbles', occurring stochastically with rate a 
[5]. The mean run length is I = v/a and duration 1/ct; 
at larger length and time scales Fick's law is obeyed, with 
diffusivity D = v 2 /da in d dimensions [5]. This model 
offers an important paradigm for a diffusion process that 
is fundamentally non-Brownian. Subtle consequences of 
this are manifest for particles in external force fields, such 
as gravity or a harmonic trap [TO]- In the first case, the 
gravitational decay length A falls strictly to zero when 
the gravitational force / exceeds the propulsive force f p , 
in contrast to Brownian particles for which A = D/f 
|10j . In a harmonic trap (/ = —kr), particles are strictly 
confined within a radius r* — f p /k; and for £<;r* the 
maximum density occurs at r ~ r* not r — 0. In this 
limit, a particle in the trap interior rapidly swims out to 
r* and stays there a long time until its next tumble [TU] • 

The qualitative physics of the aforementioned results 
is robust to both a distribution in v, or a residual true 
Brownian diffusivity. On the other hand, because there is 
no underlying free energy (which would give a Boltzmann 
distribution as the unique steady state), long-range hy- 
drodynamic interactions (HI) between the particles could 
have major consequences, even for steady-state behav- 
ior. Several computational approaches to address hydro- 
dynamics have been developed [5], but none have ad- 
dressed the basic physics problems considered below: (a) 
sedimentation in a container with a solid bottom; (b) 
confinement by a harmonic trap; and (c) Poiseuille flow 



between parallel plates. These we consider at small but 
finite particle density, so that in (a,b) only the far-field 
hydrodynamics are important. In (c), the main hydrody- 
namic effect is instead the coupling to an imposed flow. 

Problems of bacteria in force fields may appear to have 
little direct relevance to biology [5] . This could explain a 
surprising lack of experiments on both bacterial sedimen- 
tation, e.g., by centrifugation, and trapping (where the 
regime €^jr* might be achieved using recent optical meth- 
ods [H].) We argue that these simple cases demand to 
be understood before one can claim to explain more com- 
plex (and biologically relevant) ones, and hope our work 
will stimulate new experiments to help fill such gaps. 

In this Letter we address run-and-tumble systems con- 
fined by gravity, traps or walls. In the first two cases, 
our goal is to see whether the nonequilibrium steady- 
states found without HI [TU] survive with HI present. For 
sedimentation, we find only weak effects of HI. In con- 
trast, for a harmonic trap, we find that only for I <C r* 
is the near-Gaussian distribution seen without HI main- 
tained; whenever £>r* the 'density-inverted' steady state 
is destroyed by HI, replaced instead by a remarkable self- 
assembled pump-like structure. We explore in detail the 
origin of this instability, in which the local co-alignment 
of swimmers causes HI to add coherently rather than with 
random signs, vastly enhancing their effects. Thirdly, we 
address swimmers confined between walls at separations 
h. HI do not prevent swimmers from accumulating at 
the walls [7], where a weak Poiseuille flow causes strong 
upstream alignment, and hence, for h < £ a net upstream 
particle flux. This mechanism is much simpler than one 
explored previously at larger h j2]. 

An obstacle to simulations is the expense of handling 
near-field aspects of the HI [5] . These are clearly impor- 
tant at high enough density, but depend on a swimmer's 
precise geometry and stroke. In contrast, the far-field 
flow around a swimmer, in the absence of any body force 
acting on it, is of universal stresslet form [5], with force 
dipole strength s = vf p a. Here a is a hydrodynamic ra- 
dius, and v is order one; the organism is 'extensile' for 
v > 0, 'contractile' for v < 0. (Most bacteria are exten- 
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sile.) At low concentrations, modeling only the stresslet 
fields should capture any universal macroscopic conse- 
quences of HI. As shown below, this simplification dras- 
tically reduces the computational cost. 

Method: We use the lattice Boltzmann (LB) method 
to handle the fluid momentum transport [12], coupled 
to a description of each swimmer as a pair of off-lattice 
point particles at fixed separation ua (we choose v = 2). 
At these points, forces ±f p (of fixed magnitude, with + 
denoting the 'head' of the swimmer) are exerted on the 
fluid, creating an extensile stresslet. These forces are 
mapped onto the LB lattice by an optimal discretization 
|13| . alongside any external force f which we assume to 
act only on the head. (This force is passed directly onto 
the fluid: there is no particle inertia.) A local fluid veloc- 
ity (excluding a self-term) is then interpolated from the 
lattice to the head particle [13]; the swimmer co-rotates 
with the local fluid and moves relative to it with an ad- 
ditional velocity u = (f p + f)/6Trr)a. Although pointlike 
for forcing purposes, our particles have finite hydrody- 
namic radius, a — 0.05; particle scale Reynolds numbers 
are small, as required (5 x 10 -3 ). As shown in [T3] our 
method efficiently simulates dilute forced colloids, using 
a values much smaller than the lattice spacing (unity). 
Indeed with a = 0.05 we have simulated sedimentation of 
2 16 = 65536 swimmers at volume fraction 4> ~ 6.5 x 10~ D , 
on a lattice of size 128 2 x 32, using a serial code. In 
contrast, previously published works using fully resolved 
algorithms are limited to at most a few hundred parti- 
cles [5] [5] . A fully resolved LB study of the same system 
would involve a, and thus L, about 20 times larger |13j . 
This would require use of very large parallel computers. 
With the latter, we hope in future to use our far-field 
code to address situations involving millions of particles. 

Sedimentation: We now turn to our results, first briefly 
outlining those for sedimentation. (An additional figure 
and further discussion are available in appendix [Al) We 
have computed steady state density profiles p(z) under 
a constant gravity force /, at various w = f / f p , in both 
small and large systems (1000 and 65536 particles re- 
spectively), with hard walls at the base and ceiling and 
periodic boundary conditions horizontally in both cases. 
Far from a proximal region (of height ~ £ for w -C 1, 
in which the density is strongly perturbed by the basal 
wall) our numerical density profiles closely resemble the 
analytic result of |10j (exact in an unbounded domain 
without HI): p(z) ~ exp[— z/X] with X(w) —> as w 1. 
Indeed, the exponential form is maintained (within error) 
with HI present, and X(w) has for small w a very similar 
form to that without HI. However as w — > 1, X seem- 
ingly remains higher than predicted (see appendix [A]) . 
Nonetheless, the role of HI in altering the steady state re- 
mains modest: in particular, we see no evidence of HI in- 
ducing macroscopic flow patterns (which, in the absence 
of a free energy, would be possible in principle). This 
strongly contrasts with our results for traps, to which we 
turn next. 

Traps: Figjl] shows steady state densities p(r) for par- 
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FIG. 1: (Color online) Steady state density for particles con- 
fined to harmonic traps. Solid lines: LB simulations; dashed 
lines: numerics of [10] . Key: £ = £/r* values. Note that 
without HI particles cannot exceed the trap radius r*. 



tides confined to harmonic traps with various ratios of 
run length to confinement length, £ = £/r*. These are 
compared to the numerics of [TU] without HI; the latter 
were also used as the initial condition for our runs. For 
small C, we see very little effect of HI: both algorithms 
show the Gaussian distributions expected for particles of 
diffusivity D in a harmonic trap. For larger £ however, 
HI dramatically destabilize the 'density-inverted' state 
(with p[r) maximal at r*) that arises without HI. Rather 
than approaching a steady state with zero macroscopic 
flux (as the non-HI system does), the system moves to 
an attractor in which rotational symmetry is broken by 
the formation of a swarm of outward-swimming particles. 
To a good approximation the swarm remains stationary 
at a radius r s < r* , where the propulsive force balances 
the collective drag and the trapping force. The latter is 
passed on to the fluid to create a macroscopic flow: one 
has in effect a self- assembled pump. The resulting flow 
is mainly of stokeslet form, with fluid flowing opposite 
to the swimming direction. (A slow rotation of this di- 
rection, which might cancel the stokeslet term on time 
averaging, is detectable, but the rotation rate vanishes 
as a — > 0.) When a swimmer in the clump tumbles, it is 
ejected in a random direction. The velocity gradient of 
the stokeslet flow rotates it to point upstream; confined 
by the trap, it eventually must rejoin the swarm. A series 
of snapshots is given in Fig{2j and a movie in [T3] . 

This emergent structure contains regions of high den- 
sity whose details will depend on near-field physics that 
we do not resolve. However, the initial instability is well 
captured by far-field hydrodynamics. The end-state is 
likewise robust, in that any instability leading to for- 
mation of a single oriented swarm within the confines 
of a trap (harmonic or otherwise) will cause a similar 
'self-pumping' state. To better understand the initial in- 
stability, we consider the limit of infrequent tumbling, 
Ijr* ^> 1. Here, without HI, the unique steady state 
comprises a thin uniform layer of outward swimmers at 
r* , each with its propulsive force balanced by the exter- 
nal trapping force [10 . Coupling this to a solvent, we 
have locally a two dimensional layer of particles exerting 
inward point forces (stokeslets) on the fluid. 
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FIG. 2: (Color online) Simulation of 10 3 swimmers in a harmonic trap. The trap radius is indicated by the circle; £ = 8. 
Swimmers are marked with cones and shaded by local density in units of the mean (key: upper color bar). For the final image, 
randomly selected particles' trajectories over the preceding time interval At = 0.5a -1 are shown. Arrows depict fluid velocity 
colored by magnitude in units of the swim speed v (key: lower color bar). For a movie see supplementary material |14j . 




FIG. 3: (Color online) Snapshot showing the early stage in- 
stability in a simulation of 1000 swimmers, initially uniformly 
and randomly distributed on the surface of a trap. Green ar- 
rows: local flow field. Red lines: swimmer trajectories. Blue 
arrows: swimming direction. Grey arrows: radial direction. 
Swimmers in the dense patch (centre) move radially inwards. 
Those either side are initially advected away from the patch, 
but rotate their swimming direction towards it. 

Because the particles on the initial shell can have not 
only tangential displacements (causing density changes) 
but also radial and orientational ones, a formal stability 
analysis of the coupled hydrodynamics of the swimmers 
in this layer is not practicable. (Moreover, noise in the 
run and tumble dynamics could alter the conclusions.) 
Nonetheless, by a careful series of simulations (detailed, 
with additional figures, in appendixjB]) we have identified 
that the instability mode is caused initially by tangential 
density fluctuations on the surface of the trap (r = r*). 
Such density fluctuations are inevitable if the total num- 
ber of swimmers is not infinite, and lead to instability via 
the following mechanism. Any surface patch that hap- 
pens to be denser than the surrounding ones will gener- 
ate locally an excess stokeslet-like flow (not cancelled by 
the contributions from distant parts of the shell). This 
flow has two effects. First, it advects the swimmers in the 
dense patch towards the center of the trap; second, it ro- 
tates neighboring particles so that these start swimming 
towards the dense patch, creating positive feedback (Fig. 
[3]). The resulting clump creates a macroscopic flow that 
sweeps the remaining particles towards itself (Fig. [2]). 

On the other hand, if in our simulations we radially 
displace a patch without altering its density, the feed- 
back is negative; this also applies for a patch in which 
only the orientation of the swimmers is perturbed. Thus, 



although radial and orientational modes are strongly ex- 
cited in the subsequent dynamics, they are not the cause 
of the initial instability (see appendix [B]) . Notably the 
same instability is present for contractile rather than ex- 
tensile swimmers [TS]. This is fully consistent with the 
mechanism above, in which the stokeslet from the con- 
fining force is transferred to the fluid by a stationary 
swimmer (the sign of whose stresslet is then less impor- 
tant). 

Since the densest initial patch usually outstrips its 
competitors to create a dominant stokeslet, there is little 
excitation of modes above the first spherical harmonic. 
(This can also be explained by expanding in such har- 
monics the flow arising from a nonuniform shell of point 
forces [T3].) The instability reported here is somewhat re- 
lated to that found by Crowley, who showed that dense 
regions in a sedimenting suspension fall more rapidly |16j , 
with the fastest growth at long wavelengths. However, 
the initial force balance is different enough to prevent a 
direct mapping from one instability to the other. 

Upstream swimming: We turn finally to run-and- 
tumble particles confined between parallel plates at sepa- 
ration h. In the absence of HI, we expect the majority of 
particles to be near the walls whenever £>h. The reason 
is the same as for the spherical trap, but the Hi-induced 
instability of the layer found in that case is suppressed 
here by the no-slip boundary condition at the wall. 

The dynamics of a swimmer adjacent to a wall can be 
complex [2, 3, 8J. Instead of resolving the near- field HI we 
apply a truncated Lcnnard Jones (LJ) potential (range 
0.25 lattice sites) which balances the normal component 
of the propulsive force. The tangential component still 
leads to motion. In practice, near-field HI could reduce 
or enhance this motion and also can, for chiral swim- 
strokes, lead to circular orbits [3]. Treating these com- 
plexities would introduce nonuniversal parameters into 
the model; by neglecting them, we can clarify whether 
or not they are essential in the context of Poiscuillc flow. 
Here upstream swimming has been observed experimen- 
tally, and a detailed mechanism proposed that requires 
chiral swimming 2 . Without questioning this result for 
the regimes addressed in [2J, we note that for narrow 
micro-channels (£<^h), a much simpler mechanism is also 
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FIG. 4: (Color online) Main figure: density (left) and orien- 
tation (right) profiles. Parameters: u/v = 1/2,1/16,1/128 
(triangles, circles, squares) and £/h = 2, 16, 128 (green, blue, 
pink). Inset: state diagram showing contours of the average 
swimmer speed, at interval v/4. Starting from +0.5w in the 
top left (solid blue), the value decreases to zero (thick black), 
further reducing to — 0.75v on the right (dashed pink). At the 
wall, particles feel a truncated LJ potential of range 0.25. 



at work. Recall that in this regime, any swimmer in the 
bulk of the channel moves to the wall layer and stays 
there a long time before tumbling. On approach to the 
wall any such swimmer is rotated by the velocity gradient 
so as to align its swimming direction upstream. This ef- 
fect is strong whenever the product of the wall shear rate 
and the transit time across the gap is large; this equates 
to u/u<;l with u the flow speed at the mid-plane. This 
condition ensures that for a uniform p (I <C h) the up- 
stream wall flux is outweighed by a downstream bulk 
flux. However for £>h, p is peaked at the walls, and a 
net upstream flux results. For typical particle trajecto- 
ries see appendix [D] Fig [2] shows the region in parameter 
space where there is a net upstream flux, and plots of 
the particle density and mean orientation in this region. 
These results are slightly shifted quantitatively, but not 
qualitatively, by increasing the LB resolution, or adopt- 
ing slightly different (but still inelastic [TDJ) wall collision 
rules. 

Conclusion: Above we have presented results from 
an efficient hydrodynamic simulation of dilute run-and- 
tumble swimmers. Under gravity we find only perturba- 
tive effects of HI on the flux- free steady state of [ID] . 
Much stronger effects are found in geometries where, 
without HI, all particles are swimming locally in the 
same direction. This result is consistent with a crude 
power counting argument: for a sheet of swimmers, the 
1/r 2 velocity contributions from each stresslet cannot 
give a large effect unless they add coherently. This 
happens for the trap and the parallel plate geometry, 
but only when the run-length exceeds the confinement 
length. For harmonic traps, the flux-free steady state 
is then replaced by one in which a symmetry-breaking 
swarm of swimmers acts as a pump. For swimmers in 
microchannels, outward-oriented particles at the confin- 
ing walls are aligned by shear to create an upstream par- 
ticle flux without the need for chiral motion. We hope 
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FIG. 5: (Color online) Main figure: X/t vs. 1 — w from fits 
to the distal part of p(z) in a system of 1000 particles. Solid 
line: prediction of [TO]; circles, squares: simulations with HI 
at two resolutions. Inset: profiles of the swimmer density 
with the fitted exponentials (dashed). For better statistics, 
we repeated the runs for w = 0.9, 0.95, 0.975 in a 128 2 x 32 
system with 2 16 = 65536 swimmers and found extremely sim- 
ilar A values. We checked that for w > 1 the profile collapses 
entirely if we switch off HI after its formation. 



that these predictions will stimulate new quantitative ex- 
periments on the fundamental physics of self-propelled 
micro-organisms in suspension. 
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and EP/H027254 for funding. MEC is funded by the 
Royal Society. 



Appendix A: Results for sedimentation 

Fig[5]shows results for steady state density profiles p(z) 
under a constant gravity force /, at various w = f/f p , 
in a system of 1000 particles with hard walls at the base 
and ceiling and periodic boundary conditions horizon- 
tally. As discussed in the main text, far from a proximal 
region the exponential form predicted without HI in [10] 
is maintained (within error) with HI present. However as 
w — > 1, A seemingly remains higher than predicted. The 
details depend however on the resolution used (Figj5]) 
and our studies suggest that the larger fitted A is in fact 
a proximal effect caused by the random stirring of the 
fluid by particles near the wall, causing upward advec- 
tion of some of the particles above. If so, we expect 
this effect to die out at larger distances with the asymp- 
totic decay length A reverting to its non-HI value. In any 
case, the role of HI in altering the steady state appears 
modest: in particular, we see no evidence of HI inducing 
macroscopic flow patterns (which, in the absence of a free 
energy, would be possible). 
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FIG. 6: Forces on a swimmer as function of its distance to 
the center of the trap. The upper, open-headed arrows show 
the forces on the swimmer's body and the lower, filled ar- 
rows show the forces on the fluid. P is the propulsive force, 
D the drag and F ex t the external, trapping force. The very 
bottom arrows decompose the swimmer's effect on the fluid 
into stokeslet and stresslet components, illustrating that the 
external, trapping force is effectively passed onto the fluid. 



Appendix B: The origin of the trap instability 

In the following we show how the disorder in the po- 
sition of swimmers on a shell generates the instability, 
and compare it with other possible mechanisms that are 
shown to be less relevant. To do so we consider an initial 
condition where the swimmers have swum as far as they 
can and are thus located on the surface of a sphere of 
radius r* = f p /k. In the absence of any flow, the only 
forces acting on the swimmers are their propulsive forces 
and the trapping force. Because both the fluid and the 
swimmers are at rest, there is no drag force and the only 
force felt by the fluid is the propulsive force: effectively 
the external trapping force has been passed onto the fluid 
and the effect of the swimmer on the fluid is to create a 
static, inward pointing stokeslet whose amplitude equals 
the propulsive force (see Fig. |6| . 

Because our swimmers are point-like particles at zero 
Reynolds number, the flow created by N swimmers is 
simply the superposition of the flow created by each one 
of them. To understand the effect of such flows on a 
swimmer, let us consider the simpler situation of an in- 
ward pointing stokeslet of magnitude p held fixed at the 
north pole of a sphere of radius r* and a swimmer located 
at an Euler angle 0,<p — on the surface of this sphere 
(see Fig. [7|. The flow field created at this point by the 
stokeslet is given by [T7] 
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FIG. 7: We consider the flow field generated by a stokeslet p 
(red arrow) held fixed at the north pole of a sphere of radius 
r* and its effect on a swimmer located at Euler angles 9 and 
ip = 0, whose orientation (magenta arrow) makes an angle $ 
with the normal to the sphere. The vorticity is indicated in 
blue and tend to rotate the swimmer towards the stokeslet. 



For small 9, that is when the swimmer is in the vicinity 
of the stokeslet, the flow has two effects: the swimmer is 
advected inwards (u z < 0) and rotated by the vorticity 
field towards the stokeslet (u> x = uj z = while uk, < 0, 
the swimmer rotates counterclockwise on figure nh . As 
time goes on, the swimmer thus also start swimming to- 
wards the stokeslet. This argument can be made more 
quantitative by considering the dynamics of the swimmer 
(neglecting the effect of the swimmer on the flow): 
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where 7 is the effective friction coefficient of the swim- 
mer [TS], f p the propulsive force of the swimmer and 
ft the trapping force. Let us then introduce e n = 
sin 9e x + cos 9e z and eg — cos 9e x — sin 9e z the unit vec- 
tors normal and tangent to the sphere at the position 
of the swimmer, so that the position of the swimmer is 
r = ?'e n . Using the expression of the flow, the equations 
of motion of the swimmer can be rewritten 
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Let us now consider the dynamics of a swimmer that 
is initially pointing out radially ($ = 0) and is close to 
the stokeslet, say #<7r/4. Initially, one finds as expected 
r < so that the swimmers falls inwards. The second 
equation show that, initially, the flow advects the swim- 
mer away from the stokeslet (9 > at t = 0). As time 
goes on however, the angle <I> increases (third equation) 
and 9 will change sign. At larger time, the swimmer is 
indeed swimming towards the stokeslet. When the par- 
ticles move away from the surface of the trap, equations 



(B4) become more complex and we do not attempt an 
analytic approach and include these equations only to 
illustrate the qualitative behaviour of the swimmers. 
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FIG. 8: Comparison of the flow fields generated by a control 
simulation, starting from a random uniform initial condition, 
and one with the densest patch on the top. The solid lines 
correspond to (|v|oo) and are identical for the two sets of sim- 
ulations. This is expected since the initial condition with the 
densest patch on the top aims at being a simple rotation of the 
uniform one. The discrete symbols correspond to |(«)|oo- In 
this case, the control decreases much more that the simulation 
with the densest patch on the top which shows that density 
fluctuation on a scale attained by the random distribution of 
swimmers suffice to control the instability. 



Let us now come back to the problem of N swimmers 
in the absence of any external flow. Firstly, note that 
by symmetry one sees that a perfectly uniform layer of 
stokeslets cannot create any flow and is thus in mechan- 
ical equilibrium. Because we consider a finite number 
of swimmers, the density of stokeslets on the surface on 
the sphere cannot, however, be uniform but can be de- 
composed as the sum of two contributions: a perfectly 
uniform layer of inward pointing stokeslets, representing 
the average density, and a layer of inward and outward 
pointing stokeslets representing the more dense and less 
dense regions, i.e., the density fluctuations. Close to the 
denser regions, one can thus expect to see the scenario 
described above for the motion of a swimmer in an in- 
ward pointing stokeslet flow. This is indeed what we see 
in a simulation of 1000 swimmers randomly distributed 
on the boundary of the trap (see Fig. 3). Around the 
denser patch located at the top, the far-field flow created 
by the swimmers resembles that of a stokeslet. [19] Fur- 
thermore, the trajectories of the neighboring swimmers 
resemble the one described above for a single swimmer in 
a stokeslet flow. Density fluctuations within the surface 
of the trap are thus a strong candidate for the origin of 
the instability. 

To confirm this role we compare simulations where the 
densest patch of the initial condition is placed at the top 
of the trap with those obtained from random uniform ini- 
tial conditions. Let us first specify what we mean by a 
dense patch. Since there are 1000 swimmers, every patch 



of surface S contains on average N s = 250S/[n(r*) 2 ] 
swimmers, with fluctuations of the order of y/N s . We 
thus choose a patch size S that is a small fraction of the 
sphere but large enough that fluctuations are not much 
larger than the mean density. We achieve this by con- 
sidering a spherical cap defined by 6 < tt/12, which con- 
tains on average 17 swimmers. To cover the sphere with 
patches of this size, one would need roughly 60 patches, 
which means that the densest patch contains on average 
27 swimmers. [20] Our 'dense patch' initial condition is 
thus constructed as follows: each swimmer is placed uni- 
formly on the sphere with probability p or uniformly on 
the top spherical cap with probability 1 — p. In practice, 
p is chosen so that the average number of swimmers in 
the spherical cap is 27. We then want to compare the in- 
stability between this initial condition and the one where 
swimmers are uniformly randomly distributed (from now 
on, we refer to this initial configuration as 'control'). If 
tangential density fluctuations are not the cause of the 
instability, the statistical outcome of the 'dense patch' 
and 'control' simulations should be the same. 



To compare the instability, we can look at the gener- 
ated flow fields. From our lattice Boltzmann simulations 
we get a 3D grid with fluid velocities v.; at the nodes. Any 
measure of the intensity of the flow yields qualitatively 
similar results and we thus consider |v|<x, = sup i |vj|. In 
figure [8] we report (|v|<x,) and |(v)|oo, where the averages 
are done over 50 runs both for the control and the dense 
patch cases. We see that the two initial conditions give 
very similar (|v|oo): this validates our 'dense patch' ini- 
tial condition which is not very different, once rotated ap- 
propriately, from a 'control' one. Since the random uni- 
form case is perfectly isotropic, we expect the flow field to 
vanish when averaged over many runs. When comparing 
I ( v) | oo , we indeed see that such is the case for the control; 
by increasing further the number of runs over which the 
average is done, the non-zero residual flow would con- 
tinue to decrease. For the dense patch case however, the 
flow remains much larger. Because the denser patch at 
the top has to compete with other dense patches, there 
is some randomness in the orientation of the flow gener- 
ated by the instability and the intensity of the flow has 
decreased. It remains however three times larger than 
the control case, which shows that fluctuations of den- 
sity of the scale of those observed in the simulations are 
sufficient to create the instability. 



A less quantitative but more direct comparison can be 
done by starting a control simulation, locating the patch 
where the instability starts and rotating the view in such 
a way that it is at the top. One can then compare it with 
a 'dense patch' run (see figure [9]). Both configuration 
are indeed very similar, confirming the role of tangential 
density fluctuations in creating the instability. 
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FIG. 9: Comparison between a random uniform initial condition (left) and one where the densest patch is 'seeded' on the 
top (right). The random uniform one has been rotated so that the instability starts on the top. We see very similar flow 
patterns and displacement of swimmers close to the north pole. Note that because of secondary dense patches, the rest of the 
configuration can slightly differ. 
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FIG. 10: Comparison between control, dense patch, and ra- 
dially displaced initial conditions for displacements ranging 
from 10% to 50% of the trap radius. Solid lines correspond 
to (|v|oo) whereas symbols correspond to |(v)|oo. All initial 
conditions produce similar flow (solid line), and thus similar 
instability. The averaged flows of the displaced initial con- 
ditions and the control are the same at large times, showing 
that these perturbations have a much smaller effect on the 
setup than the one due to density fluctuations. Note that the 
non-zero initial flow is caused by particles swimming back to 
the border of the trap. 



order'), and second they may not all be strictly pointing 
outwards ('angular disorder'). We now show below that 
these fluctuations are less relevant than the tangential 
density fluctuations. 



a. Radial Displacement 

To study the role of radial displacement, we place the 
swimmers at random over the whole sphere uniformly 
and then move downwards those in the same spherical 
cap as previously contained the dense set of bacteria. 
(There is no longer an excess density in this cap, how- 
ever.) Again we compare the instability by looking at 

' " <|v| 
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the flow fields. As can be seen in figure 
left unchanged for displacements up to half the radius of 
the trap. This means that displacements of swimmers, 
even much larger those spontaneously observed in simu- 
lations, do not seem to generate any kind of instability 
with a stronger effect than the one seen in their absence. 
Let us now turn to the comparision of |(v)|oo to look 
for systematic effects. At large times, the flow field is 
the same as the control one and much smaller than that 
of the dense patch; this perturbation does not fix the 
direction of the instability. Hence radial disorder, even 
for large displacements, does not control the instability 
observed from the random initial configurations. 



1. Other type of perturbations 



b. Rotation 



In addition to the non-uniformity of the density of 
swimmers on the surface of the trap, other sources of per- 
turbation could be responsible for the instability. First, 
the swimmers can fluctuate in radial position ('radial dis- 



To investigate the effect of disorder in the swimmers 
orientation, we first place the swimmers at random over 
the whole surface uniformly and then modify the orien- 
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FIG. 11: Comparison between control, dense patch and ran- 
domly rotated initial conditions for several rotation angles. 
Solid lines correspond to (|v|oo) whereas symbols correspond 
to |(v)|oo. All initial conditions produces similar strength of 
flow, and thus of instability. The average flow of the random 
rotations equals that of the control, showing that these per- 
turbation have a much smaller effect on the setup than the 
density fluctuations. 



FIG. 12: Comparison between control, dense patch, and uni- 
formly rotated initial conditions for several rotation angles. 
Solid lines correspond to (|v|oo) whereas symbols correspond 
to | (v) | oo . All initial conditions produce similar strength of 
flow, and thus of instability. The average flow of the uniformly 
rotated initial condition is however of the same order as the 
control one, showing that these perturbation have a much 
smaller effect on the setup than the density fluctuations. 



tation of the swimmers that are in the top spherical cap 
defined by < tt/12. We either generate randomly a 
rotation matrix and apply it to all these simmwers (uni- 
form rotation) or we pick up one rotation matrix per 
swimmer (random rotation) . The strength of the pertur- 
bation is controlled by varying the average rotation angle 
of the distribution of the rotation matrix. In both cases, 
(l v |oo) perfectly overlaps with the control (see figure 



11 



and 12 1, while the quantity |(v)|oo is as small as the con- 
trol. These results confirm that the initial rotation has 
no part to play in triggering the instability. 




Appendix C: Brownian Motion 

Interestingly, our stability analysis allows us to eval- 
uate the impact of Brownian motion on the trap in- 
stability. For E. coli, thermal diffusivity is D ~ 
0.3 /im 2 /s whereas its rotational counterpart is D rot ~ 
0.15 rad 2 /s [T5] . Over the time scale of a second, which is 
the average time between two tumbles, this corresponds 
to diffusion length and angle of about half a micrometer 
and 20 degrees. The simulations presented in the previ- 
ous subsections shows that these perturbations are negli- 
gible when compared to the density fluctuations. Brow- 
nian motion should thus not affect the trap instability. 



FIG. 13: (Color online) 40 randomly chosen swimmers (of 
1000 simulated) in a Poiseuille flow between plates. Swimmers 
are marked with cones and shaded by relative local density 
(key: upper color bar). Trajectories over the preceding time 
interval At = 0.5a -1 are shown, colored by magnitude of 
lab- frame speed in units of v (key: lower color bar). The flow 
profile is shown on the left (same key). 



Appendix D: Channel Flow 

An additional figure for channel flow, showing swim- 
mer trajectories, is given as Fig. |13| 
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